// Magic Software, Inc.
// http://www.magic-software.com
// Copyright (c) 2000, All Rights Reserved
//
// Source code from Magic Software is supplied under the terms of a license
// agreement and may not be copied or disclosed except in accordance with the
// terms of that agreement.  The various license agreements may be found at
// the Magic Software web site.  This file is subject to the license
//
// FREE SOURCE CODE
// http://www.magic-software.com/License/free.pdf

#include "MgcBisect1.h"

//---------------------------------------------------------------------------
MgcBisect1::MgcBisect1 (Function oF, int iMaxLevel, MgcReal fTolerance)
{
    m_oF = oF;
    m_iMaxLevel = iMaxLevel;
    m_fTolerance = fTolerance;
}
//---------------------------------------------------------------------------
bool MgcBisect1::Bisect (MgcReal fX0, MgcReal fX1, MgcReal& rfRoot)
{
    // test two endpoints
    MgcReal fF0 = m_oF(fX0);
    if ( MgcMath::Abs(fF0) <= m_fTolerance )
    {
        rfRoot = fX0;
        return true;
    }

    MgcReal fF1 = m_oF(fX1);
    if ( MgcMath::Abs(fF1) <= m_fTolerance )
    {
        rfRoot = fX1;
        return true;
    }

    if ( fF0*fF1 > 0.0 )
        return false;

    for (int iLevel = 0; iLevel < m_iMaxLevel; iLevel++)
    {
        MgcReal fXm = 0.5*(fX0+fX1);
        MgcReal fFm = m_oF(fXm);
        if ( MgcMath::Abs(fFm) <= m_fTolerance )
        {
            rfRoot = fXm;
            return true;
        }

        if ( fF0*fFm < 0.0 )
        {
            fX1 = fXm;
            fF1 = fFm;
        }
        else if ( fF1*fFm < 0.0 )
        {
            fX0 = fXm;
            fF0 = fFm;
        }
    }

    return false;
}
//---------------------------------------------------------------------------

#ifdef BISECT1_TEST

#include "MgcRTLib.h"

// Bisection should produce a root at (r2)
// Output:
// root (+0.707031)
// func (-0.000076)

static const MgcReal gs_fRootHalf = MgcMath::Sqrt(0.5);

MgcReal F (MgcReal fX) { return fX - gs_fRootHalf; }

int main ()
{
    MgcBisect1 kBisector(F,16,1e-04);
    MgcReal fRoot;
    if ( kBisector.Bisect(0.0,1.0,fRoot) )
    {
        cout << "root at x = " << fRoot << endl;
        cout << "function F(x) = " << F(fRoot) << endl;
    }
    else
    {
        cout << "root not found" << endl;
    }

    return 0;
}

#endif
